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1  Accomplishments/New  Findings 

•  We  have  extended  the  subpixel  smoothing  Finite-Difference  Time-Domain  (FDTD)  method  to 
material  interface  between  dielectric  and  dispersive  media  by  local  coordinate  rotation.  Our 
method  is  equivalent  to  the  previously  proposed  subpixel  smoothing  method  for  dielectric 
interface,  and  the  extension  to  dispersive/dielectric  interface  does  not  require  split  fields  so 
our  method  has  improved  the  efficiency  in  comparison  to  the  previous  proposed  split  held 
approach. 

•  A  novel  stable  anisotropic  FDTD  algorithm  based  on  the  overlapping  cells  has  been  devel¬ 
oped  for  solving  Maxwell’s  equations  of  electrodynamics  in  anisotropic  media  with  interface 
between  anisotropic  dielectrics  and  dispersive  medium  or  Perfect  Electric  Conductor  (PEC). 
The  previous  proposed  conventional  anisotropic  FDTD  methods  suffer  from  the  late-time  in¬ 
stability  due  to  the  extrapolation  of  the  held  components  near  the  material  interface.  Our 
anisotropic  Overlapping  Yee  (OY)  FDTD  method  is  stable,  as  it  relies  on  the  overlapping  cells 
to  provide  the  collocated  held  values  without  any  interpolation  or  extrapolation.  Numerical 
results  on  eigenvalue  analysis  confirm  that  our  method  is  stable.  We  has  applied  our  method 
to  simulate  the  electromagnetic  invisibility  cloaking  devices. 

•  We  have  extended  our  recently  proposed  OY  FDTD  method  to  locally  non-orthogonal  grids, 
with  application  to  the  optical  force  computation  on  nanoparticles.  The  subpixel-smoothing 
FDTD  technique  has  successfully  achieves  second-order  accuracy  by  using  an  inverse  dielectric 
tensor.  However,  a  remaining  challenge  is  to  accurately  handle  objects  with  sharp  corners, 
where  the  accuracy  is  still  less  than  second-order.  Our  approach  has  attained  second-order 
convergence  when  sharp  corners  present. 

•  We  have  developed  a  moving  window  full  Maxwell  solver  algorithm  with  perfectly  matched 
absorbing  layer  (PML)  boundary  conditions  in  order  to  accurately  simulate  the  propagation 
of  localized  waves  over  a  very  long  distance  (millions  of  wavelength)  in  complex  media.  An 
existing  finite  difference  moving  frame  method  developed  more  than  a  decade  ago  is  inade¬ 
quate  due  to  low  order  transparent  boundary  conditions.  Our  method  enables  the  realistic 
and  predictive  simulations  of  high  intensity  optical  pulses  in  regime  for  which  current  direct 
Maxwell  solvers  are  inapplicable  due  to  memory  and  CPU  requirements. 

•  We  have  implemented  the  Adaptive  Mesh  Refinement  (AMR)  FDTD  method  to  study  the 
ground  penetrating  radar  devices.  This  work  was  included  in  the  MS  thesis  of  a  graduate 
student. 


2  Subpixel  smoothing  FDTD  Maxwell  solver  for  dielectric/dispersive 
interface 

We  have  developed  a  new  way  to  extend  the  subpixel  smoothing  FDTD  algorithm  proposed  in  [1,  2] 
to  material  interfaces  between  dielectric  and  dispersive  media  by  using  local  coordinate  rotation. 

The  advantage  of  our  method  is  the  efficiency  as  it  does  not  require  the  storage  and  computation 
of  the  split  fields  and  additional  equations. 

The  constitutive  equation  that  converts  the  electric  flux  D  to  electric  held  E  is  given  by 

E  =  €o  1r1D,  (1) 


2 


Figure  1:  Electric  fields  near  the  material  interface.  Ex  and  Ey  are  the  electric  fields  in  Cartesian 
coordinates.  Ejy  and  Et  correspond  to  the  electric  fields  perpendicular  and  tangential  to  the 
material  interface  respectively. 


where  eo  is  the  permittivity  in  vacuum  and  e-1  is  the  inverse  permittivity  tensor.  Our  coordinate 
rotation  algorithm  solve  it  in  the  following  steps: 

First,  we  compute  the  normal  and  tangential  components  of  the  electric  fluxes  Dn  and  Dt  using 
coordinate  rotation  (Fig.  1) 


where  R  is  the  rotation  matrix 


(2) 

(3) 


Second,  we  compute  the  normal  and  tangential  components  of  the  electric  fields  En  and  Et  by 
solving  the  constitutive  equation  in  the  rotated  coordinates. 


(  En  \  =  1_  f  <  e-1  >  0  \  /  Dn  \ 

V  Et  j  eo  V  0  <  e  >-x  )  V  Dt  J  ’ 


(4) 


where  <  •  >  denotes  the  volume  average.  In  the  rotated  coordinates,  the  constitutive  equation  is  a 
diagonal  system  so  it  is  solved  as  two  separated  equations: 


Ejv 


—  <  e  1  >  Djv, 
eo 


(5) 


and 


Et  =  —  <  e  >  1  Dt , 
eo 


(6) 


Finally,  we  update  the  elelctric  fields  in  the  regular  Cartesian  coordinates  Ex  and  Ey  using 
inverse  rotation 


Ev 


=  R 


-l 


En 

Et 


(7) 


Combing  three  steps,  the  constitutive  equation  becomes 


Ex 

En 


=  —R 
eo 


-l 


<  e  1  > 


0 

<  e  > 


-i 


R 


Dx 

Dn 


(8) 
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Let  P  =  NtN  be  the  projection  matrix  where  N  is  the  unit  normal  vector,  we  see  that  equation 
(8)  is  equivalent  to  equation  (1)  if  we  let 

r1  =  P  <  e”1  >  +(1  -  P)  <  e  >~1,  (9) 


where  e"1  is  the  inverse  permittivity  tensor  for  the  subpixel  smoothing  algorithm  [1],  Therefore,  we 
see  that  the  coordinate  rotation  results  in  the  same  formulation  as  the  subpixel  smoothing  method 
for  a  dielectric  interface. 

To  illustrate  how  to  solve  equations  (5)  and  (6)  on  the  material  interface  between  dispersive 
and  dielectric  media  without  using  split  fields,  we  use  the  Lorentz  dispersive  model  as  example.  As 
shown  in  Fig.  1,  we  assume  that  ei  is  a  constant  (dielectric  material)  and  e2  represents  a  Lorentz 
dispersive  medium  with  dielectric  function  given  by 


62  (w) 


^OO 


at 

< jo 2  —  i'yuj  —  u2 


(10) 


In  equation  (6),  we  first  compute  volume  averaged  permittivity: 


<  e  >=  (3e2{u)  +  (1  -  (3)e i, 


(11) 


where  (3  represents  the  ratio  of  the  shaded  area  of  the  Fig.  1  to  the  whole  area  of  the  figure.  After 
some  simple  calculation,  we  get 


aoj~ 


<  e  >=  eoo  - 


OO  o  .  ^  ^  9  ) 


(12) 


where  =  /3eoo  +  (1  —  (3)e i,  a  =  (3a ,  Cjp  =  lop ,  and  7  =  7.  By  introducing  the  polarization  Pt,  we 
have 

Dt  =  eo  <  e  >  Et  =  eo(eoo Et  +  Pt),  (13) 


where 


Pt  =  ~ 


aut, 


uj2  —  i^foj  —  uj2 


Et- 


(14) 


Equation  (14)  is  a  Lorentz  model  and  can  be  solved  by  the  Auxiliary  Differential  Equation  (ADE) 
algorithm  [3]: 

dfr2  +  7^  +  tf,pT  =  uC?vEt.  (15) 

For  the  normal  component  of  the  electric  field  in  equation  (5),  we  have  the  harmonic  averaging 


<e  1  >—  (3 — rr  +  (l-^)  — , 
e2(w)  ei 


which  leads  to 


<  e”1  > 


aujp 

ujz  —  xoj'y  — 


(16) 


(17) 


,2  /c2  “2  _  ,  ,2  1  ,.,2( 

“1  /  S  7  ^ p  + 

the  previous  case,  from  equation  (5),  we  get 


where  Coo  =  CooCi/^,  a  =  a(3e2/£2,  u2  =  uj2+  u2(l  -  (3)a/£ ,  and  £  =  /3ei  +  (1  -  (3) e qq.  Similarly  to 


L>at  =  eo  wry  — ffiv  =  60(^00 -E)v  +  Dv), 


(18) 
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(a)  (b) 

Figure  2:  (a)  Electric  fields  intensity  along  the  x-axis  passing  through  the  center  of  a  cylinder,  (b) 
Relative  error  versus  resolution. 

where  P/v  is  the  normal  component  of  the  polarization 

aui 

PN  = - n - “ - 3“o  P/v 

U+  —  IJUJ  —  Ulp 

Equation  (19)  is  also  a  Lorentz  model  and  can  be  solved  in  the  same  way  as  the  previous  case  in 
equation  (14). 

The  split  field  subpixel  method  proposed  in  [4]  consider  a  general  interface  between  two  disper¬ 
sive  media.  If  one  side  of  the  interface  is  a  dielectric  medium,  it  can  be  reduced  to  two  auxiliary 
variables  since  <  e-1  >  in  equation  (9)  is  given  in  equation  (17).  However,  split  fields  are  still 
required  since  it  cannot  be  further  simplified,  unless  by  using  a  common  denominator  which  will  re¬ 
sult  in  higher  order  differential  equations  due  to  the  higher  order  u>  terms.  In  contrast,  our  method 
is  a  diagonal  system  so  the  constitutive  equation  can  be  solved  separately  and  the  field  splitting  is 
avoided.  A  limitation  of  our  method  is  that  it  applies  to  the  case  where  one  side  of  the  interface  is 
a  dispersive  medium.  For  an  interface  between  two  dispersive  media,  split  fields  are  still  required. 
We  would  like  to  point  out  that  the  normal  and  tangential  components  of  the  electric  fields  and 
fluxes  ( En ,  Et,  Thv,  and  Dt)  in  our  algorithm  serve  as  intermediate  variables,  so  except  for  a  few 
temporary  variables,  no  additional  computational  memory  is  required. 

We  test  our  method  (the  coordinate  rotation  (CR)  subpixel  smoothing)  for  a  scattering  problem 
and  compare  it  with  the  split  held  (SF)  subpixel  smoothing  method,  the  standard  FDTD  method, 
and  the  Mie  solution.  In  our  numerical  simulation,  a  cylindrical  particle  with  refractive  index 
n  =  0.23  + 2. 97i  is  modeled  by  the  Lorentz  model.  Fig.  2  shows  the  electric  held  intensity  along  the 
x-axis  passing  through  the  center  of  the  cylinder  and  the  comparison  of  the  relative  error  versus 
the  resolution.  The  subpixel  smoothing  method  consistently  achieves  smaller  error  in  comparison 
to  the  staircased  FDTD  results.  We  also  see  that  although  our  method  does  not  achieve  second 
order  accuracy,  it  does  signihcantly  improve  the  accuracy.  Table  1  shows  the  comparison  on  CPU 
time  and  memory  usage  for  the  three  methods.  In  comparison  to  the  SF-subpixel  smoothing,  our 
method  appears  to  have  better  accuracy  when  the  numerical  grid  resolution  is  low  and  it  saves 
about  10%  of  memory  and  about  25%  of  CPU  time. 


(19) 
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Table  1:  Comparison  of  CPU  time  and  memory  usage 


Method 

A  = 
Memory 

A/40 

Runtime 

A  = 
Memory 

A/80 

Runtime 

staircase 

13  MB 

7  sec 

50  Mb 

77  sec 

CR 

21  MB 

9  sec 

85  Mb 

120  sec 

SF 

23  MB 

12  sec 

93  Mb 

150  sec 

3  Overlapping  Yee  FDTD  method  for  Material  Interfaces  between 
Anisotropic  Dielectrics  and  General  Dispersive  or  PEC  Media 


The  anisotropic  FDTD  method  have  been  studied  for  many  years  [5,  6,  7,  8,  9,  10,  11,  2,  12]. 
Reference  [11]  pointed  out  that  for  stable  algorithm  it  is  sufficient  to  have  the  finite-difference 
operator  that  converts  D  to  E  (the  material  matrix)  to  be  symmetric  and  positive  semidefinite. 
In  the  same  paper,  the  authors  developed  a  stable  FDTD  algorithm  for  anisotropic  dielectrics  and 
showed  that  the  previously  developed  methods  lead  to  asymmetric  material  matrices,  which  implied 
instability.  Stability  analysis  in  case  of  non-uniform  dielectric  was  provided.  The  method  presented 
in  [11]  only  consider  the  case  where  the  anisotropic  dielectrics  are  away  from  the  material  interface. 
When  the  anisotropic  dielectrics  are  very  close  to  the  material  interface,  such  as  the  interface 
between  dielectrics  and  dispersive  media  and  between  dielectrics  and  Perfect  Electric  Conductor 
(PEC)  boundaries,  extrapolation  is  required  [7]. 

For  anisotropic  materials,  the  constitutive  equations  are  given  by 


Bx  \  /  t^xx  l^xy  /  tr z  \  /  Hx  \ 

By  J  —  /Xo  I  f^yx  M1 yy  l^yz  J  I  By  J  ,  (21) 

Bz  /  \  f^zx  l^zy  l-^zz  /  \  Hz  ) 

where  matrices  e  =  {e^}  and  /r  =  are  generally  non-diagonal.  In  order  to  solve  equation  (20) 
for  the  electric  fields  E,  all  components  of  the  electric  fields  and  fluxes  need  to  be  collocated  at  the 
same  location.  Similarly,  for  equation  (21),  all  components  of  the  magnetic  fields  must  be  collocated. 
The  standard  FDTD  Yee  lattice  is  staggered,  so  only  one  component  of  the  field  is  provided  at  each 
location.  In  order  to  obtain  the  other  two  components,  a  conventional  way  is  to  average  neighbor 
values.  However,  as  it  has  been  pointed  out  in  [11],  direct  averaging  causes  late-time  instability 
and  a  stable  interpolation  algorithm  has  been  proposed  in  the  same  article.  A  limitation  of  the 
method  proposed  in  [11]  is  that  it  considers  only  the  case  where  the  anisotropic  dielectrics  are  away 
from  the  material  interface.  If  the  anisotropic  dielectrics  are  very  close  to  the  material  interface, 
such  as  the  interface  between  dielectrics  and  dispersive  material  or  PEC  boundaries,  extrapolation 
is  required.  We  have  found  that  the  extrapolation  will  result  in  late-time  instability,  similar  to  the 
previous  case  where  no  material  interface  presents.  An  alternative  anistropic  FDTD  using  Finite 
Elements  near  the  PEC  boundaries  is  proposed  in  [12],  but  it  is  limited  to  PEC  boundaries  and 
does  not  handle  problems  with  magnetic  anisotropy. 

In  order  to  overcome  this  late-time  instability,  we  propose  a  new  stable  FDTD  Maxwell  solver 
for  anisotropic  media,  based  on  the  idea  of  using  overlapping  cells.  The  usage  of  overlapping  cells 
guarantees  the  collocation  of  field  components  so  that  the  interpolation  and  extrapolation  that 
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Figure  3:  A  3D  computational  cell  and  its  corresponding  overlapping  cell  in  x-y  plane. 


cause  instability  is  avoided.  This  material  interface  could  be  the  interface  between  dielectrics  and 
dispersive  medium,  the  interface  between  dielectric  and  PEC  boundary  or  possibly  other  types  of 
boundaries.  Since  the  OY  FDTD  has  the  same  numerical  discretization  as  the  standard  Yee  FDTD 
method  (but  on  multiple  grids),  the  same  source  conditions  and  perfectly  matched  layer  (PML) 
boundary  conditions  apply. 

As  shown  in  Fig.  3,  an  overlapped  mesh  is  constructed  by  shifting  the  primary  grid  by  half  cells 
in  the  x-y  plane.  The  z-coordinates  do  not  change.  Note  that  in  3D,  the  overlapped  mesh  is  not  the 
dual  mesh  since  the  dual  mesh  is  generated  by  shifting  the  grid  by  half  cell  in  all  three  directions. 
Similarly,  on  the  y-z  and  z-x  planes,  we  construct  two  other  overlapped  meshes.  Therefore,  the  3D 
OY  mesh  contains  four  Yee  lattices. 

On  each  Yee  grid,  the  standard  Yee  FDTD  method  is  applied  to  update  electric/magnetic  fluxes 
using  the  information  of  the  nearby  magnetic/electric  fields.  After  updating  the  fluxes,  the  local 
constitutive  equations  are  used  to  solve  for  electric/magnetic  fields.  For  instance,  to  update  E  at 
( i  +  2 , j,  k)  at  time  step  n,  we  solve  the  following  3x3  linear  system: 

/  D™(i  +  \,j,k)  \  /  exx  exy  exz  \  (  E^(i  +  ^,j,k)  \ 

Dy(i  +  2, J,  k)  I  =  e0  {  eyx  eyy  eyz  E™(i  +  \,j,k)  ,  (22) 

V  D™(i  +  \,j,k)  )  \  ezx  ezy  ezz  j  \  E%{i  +  \,j,k)  ) 

where  D” (?,  +  |,  j,  k),  Dy(i  +  j,  k ),  and  D™(i  +  |,  j,  k )  belong  to  three  different  Yee  lattices. 

We  have  applied  our  method  to  simulate  the  cylindrical  cloak  proposed  in  [13].  In  the  cloaking 
shell,  the  electric  permittivity  and  magnetic  permeability  are  functions  of  the  radius  r: 


er(r) 

te(r) 


Pz(r) 


r  —  R\ 

5 

r 

r 

r-Ri 

r-Rif  R-2  V 
r  \R2-Ri)  ’ 


(23) 

(24) 

(25) 


where  R2  and  R\  are  the  radii  of  the  outer  and  inner  circles  of  the  cloaking  shell  respectively,  as 
shown  in  Fig.  4.  The  inner  circle  is  a  PEC  shell.  We  consider  the  transverse-electric  (TE)  polarized 
electromagnetic  wave  propagation  in  x  direction. 

In  our  numerical  simulations,  the  incident  wave  is  a  plane  wave  propagating  in  x-direction  with 
wavelength  A  =  0.15  m.  The  object  to  be  cloaked  is  a  circular  cylinder  with  radius  R±  =  0.1  m  and 
the  cloaking  shell  has  thickness  0.1  m,  so  R2  =  0.2  m.  The  computational  domain  is  surrounded  by 
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PML 


Figure  4:  Computational  domain  for  the  FDTD  metamaterial  cloaking  simulation. 


PML  absorbing  boundaries.  To  avoid  the  singularity  at  r  =  R\,  the  cloaking  shell  is  approximated 
as  a  ring  with  eight  discrete  layers  of  homogeneous  dielectric  parameters,  similar  to  [13]  (Fig.  4). 

We  have  investigated  the  stability  of  the  conventional  anisotropic  FDTD  and  our  OY  method 
by  computing  the  eigenvalues  of  the  fully  discrete  problem  in  two  dimensions.  As  shown  in  Fig.  5, 
the  conventional  anisotropic  FDTD  method  generates  a  few  eigenvalues  outside  the  unit  circle.  As 
a  result,  the  solution  grows  slowly  and  eventually  blows  up.  In  contrast,  the  OY  method  has  all 
eigenvalues  located  on  the  unit  circle  which  guarantees  the  stability. 

Using  the  cloaking  simulation,  we  compare  the  conventional  anisotropic  FDTD  and  our  anisotropic 
OY  FDTD  methods.  Fig.  6  shows  our  numerical  simulation  results  on  the  electric  field  distribution 
near  the  cloaked  objects  after  a  few  thousands  of  steps.  Both  methods  give  correct  answer  at 
early  stage,  but  the  conventional  FDTD  method  results  in  some  spikes  near  the  inner  circle  of  the 
object,  which  is  due  to  the  extrapolation  and  leads  to  the  late-time  instability.  These  spikes  are 
amplified  as  time  steps  increase.  Our  simulations  show  that  usually  after  a  few  tens  of  thousands 
steps  (sometimes  just  a  few  thousands  of  steps),  the  solution  blows  up.  In  contrast,  the  new  OY 
algorithm  gives  smooth  and  stable  solution  for  long  time  simulation. 

4  Locally  non-orthogonal  OY  FDTD  method  and  optical  force 
computation 

The  subpixel-smoothing  FDTD  technique  proposed  in  [1]  achieves  second-order  accuracy  by  using 
an  inverse  dielectric  tensor,  and  this  method  has  been  extended  to  anisotropic  media  [2],  combined 
with  a  recently  proposed  stable  FDTD  scheme  in  anisotropic  media  [11].  However,  a  remaining 
challenge  is  to  accurately  handle  objects  with  sharp  corners,  where  the  accuracy  is  still  less  than 
second-order  [1] .  Another  way  of  eliminating  staircasing  is  to  extend  the  FDTD  methods  to  globally 
non-orthogonal  meshes  [14,  15,  16,  17],  but  these  methods  suffer  from  the  late-time  instability 
due  to  the  non-positive  definite  property  of  the  material  (projection)  matrices  introduced  by  the 
numerical  methods.  Our  recently  developed  non-orthogonal  overlapping  Yee  (OY)  FDTD  method 
has  successfully  overcomed  the  late-time  instability  problem  [18]. 

We  have  further  investigated  the  non-orthogonal  OY  method  to  model  the  optical  force  compu¬ 
tation.  Second  order  convergence  was  achieved  even  when  sharp  corner  presents.  The  OY  method 
requires  multiple  Yee  grids  which  lead  to  more  memory  and  CPU  costs.  To  improve  the  efficiency, 


(a)  (b) 


(c)  (d) 

Figure  5:  Distribution  of  the  eigenvalues  on  the  complex  plane  for  cloaking  simulation  using  the 
conventional  FDTD  method  (upper  row)  and  OY  FDTD  method  (lower  row).  On  the  right  are  the 
enlarged  figures. 


Figure  6:  Numerical  results  on  electromagnetic  metamaterial  cloaking  simulations.  Left  column: 
the  conventional  anisotropic  FDTD  method;  Right  column:  the  anisotropic  Overlapping  Yee  FDTD 
method.  Upper  row:  contour  plots  of  the  electric  field  Ey  distribution;  Lower  row:  surface  plots  of 
the  electric  field  Ey  distribution. 
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Figure  7:  Non-orthogonal  quadrilateral  meshes  near  a  tilted-square  (Left)  and  triangular  wedge 
(Right).  The  sides  of  the  particles  (in  red)  go  through  the  diagonal  vertices  of  the  grid  cells. 


we  propose  a  locally  non-orthogonal  OY  technique.  In  our  new  implementation,  the  computational 
cells  are  non-orthogonal  and  overlapped  only  in  a  small  region  near  the  curved  geometry,  and  the 
rest  of  the  computational  domain  is  regular  Yee  grid. 

We  use  the  Maxwell  stress  tensor  formulation  to  compute  the  total  optical  force  acting  on  a 
particle  due  to  a  time-harmonic  electromagnetic  field: 

<  F  >=  JJ  <  T  >  •  dS,  (26) 

where  T  is  the  stress  tensor 

Tij  =  eEiEj  +  fiH{E[j  —  —  (eE2  +  nH2)5ij.  (27) 

The  angle  brackets  indicate  time-averaged  values.  S'  is  a  surface  that  enclose  the  object.  A  simple 
choice  of  the  surface  S  is  a  cubic  box  outside  the  object. 

The  diagonal  split-cell  model  is  employed  to  construct  quadrilateral  meshes  that  avoids  the 
permittivity  averaging.  The  diagonal  split-cell  model  has  been  discussed  in  Taflove’s  book  [3]  for 
orthogonal  Cartesian  grid.  We  have  extended  this  model  to  non-orthogonal  grids.  As  shown  in 
Fig.  7,  the  quadrilateral  mesh  is  constructed  in  such  a  way  that  the  material  interface  does  not 
go  along  cell  edges  but  passes  through  the  diagonal  vertices  of  the  cells  using  a  smooth  circle 
mapping  method  similar  to  the  method  proposed  in  [19].  Besides  the  simple  test  cases  (e.g.,  a 
circular  cylinder  or  a  rectangle),  our  method  can  be  applied  to  more  complicated  geometries,  such 
as  a  triangular  wedge,  as  shown  in  Fig.  7.  By  placing  the  magnetic  fields  at  the  cell  centers  and 
the  electric  fields  along  cell  edges,  this  approach  guarantees  that  no  line  integral  of  the  electric 
field  crosses  the  material  interface  so  that  the  permittivity  averaging  is  avoided.  This  approach  is 
identified  as  the  Split-Cell  Overlapping  Yee  (SC-OY)  method.  The  split-cell  mesh  construction  and 
some  preliminary  result  of  computing  force  on  cylindrical  particle  were  presented  in  our  previous 
published  work  [20].  The  non-orthogonal  quadrilateral  split-cell  mesh  for  structures  with  sharp 
corners,  such  as  a  titled-square,  is  shown  in  Fig.  7.  The  SC-OY  method  has  at  least  two  advantages: 
(1)  it  avoids  permittivity  averaging,  so  that  the  implementation  is  simplified;  (2)  it  avoids  tiny 
and  near  180  degree  angles  for  structures  with  large  curvature  so  that  the  local  error  is  smaller. 
A  limitation  of  the  split-cell  approach  is  that  it  assumes  nonmagnetic  medium  (fj,  is  constant 
throughout  the  computational  domain). 
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Figure  8:  Locally  non-orthogonal  quadrilateral  mesh  near  a  tilted-square. 


In  general,  the  2D  OY  technique  requires  two  set  of  Yee  grids  to  be  overlapped  to  each  other,  so 
it  doubles  the  computational  cost.  In  3D,  it  quadruples  the  cost.  To  improve  the  efficiency  of  OY 
method,  we  have  employed  a  locally  non-orthogonal  approach.  That  is,  in  our  implementation,  the 
computational  cells  are  only  locally  non-orthogonal  in  the  region  that  near  the  curved  geometry. 
The  rest  of  the  domain  is  orthogonal  rectangular  grids.  A  sample  mesh  is  shown  in  Fig.  8.  In 
the  non-orthogonal  region,  the  OY  method  is  applied  and  the  standard  FDTD  method  is  applied 
elsewhere.  As  shown  in  this  figure,  the  non-orthogonal  mesh  is  only  about  10%  of  the  whole 
computational  domain,  so  that  the  overall  computational  cost  decreases  from  200%  to  110%.  In  3D, 
this  locally  non-orthogonal  technique  will  further  improve  the  performance  of  numerical  simulation. 

In  our  numerical  examples,  we  have  applied  the  nonorthogonal  OY  technique  to  compute  the 
optical  force  on  a  tilted-square  with  sharp  corner  of  90  degrees.  The  side  length  of  the  square 
is  240  nm  and  it  is  titled  for  30  degrees.  The  incident  plane  wave  has  wavelength  A  =  600  nm. 
and  propagates  in  x-direction.  The  computational  domain  is  surrounded  by  the  uniaxial  perfectly 
matched  layer  (UPML)  boundaries  in  all  directions. 

Fig.  9  shows  the  relative  errors  of  the  computed  optical  forces  on  dielectric  and  metallic  particles. 
The  numerical  result  at  very  fine  mesh  (A  =  A/400)  is  used  as  the  exact  solution.  As  shown  in 
Fig.  9,  for  both  dielectric  and  metallic  particles,  the  FDTD  method  converges  linearly,  while  the 
OY  method  is  second-order  accuracy.  For  grid  size  A  =  A/200,  the  relative  errors  of  the  OY  results 
are  about  one  order  of  magnitude  smaller  than  the  FDTD  results.  Our  numerical  results  also  show 
that  the  OY  solution  with  A  =  A/80  and  the  FDTD  solution  with  A  =  A/200  have  comparable 
errors  (about  0.5%),  but  the  OY  method  requires  only  half  resource  in  memory  usage  and  less  than 
25%  resource  in  CPU  usage. 

5  Moving  Window  FDTD  with  PML  boundaries 

We  have  developed  a  novel  moving  window  full  Maxwell  solver  algorithm  with  perfectly  matched 
absorbing  layer  (PML)  boundary  conditions  in  order  to  accuratelly  simulate  the  propagation  of 
localized  waves  over  a  very  long  distance  (millions  of  wavelength)  in  complex  media.  Our  method 
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(a) 


(b) 


Figure  9:  Relative  error  of  forces  versus  resolution  for  tilted-square  illuminated  by  TEZ  plane  wave, 
(a)  dielectric  medium  (t2  =  9);  (b)  dispersive  medium. 


enables  the  realistic  and  predictive  simulations  of  high  intensity  optical  pulses  in  regime  for  which 
current  direct  Maxwell  solvers  are  inapplicable  due  to  memory  and  CPU  requirements.  An  existing 
finite  difference  moving  frame  method  developed  more  than  a  decade  ago  is  inadequate  due  to  low 
order  transparent  boundary  conditions.  At  present,  the  wave  propagation  in  a  long-distance  regime 
are  modeled  via  asymptotically  reduced  unidirectional  models  relying  on  the  assumption  that  all 
electromagnetic  propagation  modes  are  known  and  decomposition  into  forward-backward  moving 
waves  is  possible.  Both  of  these  assumptions  fail  in  the  high  contrast  and  nonlinear  media,  and  for 
ultra-short  pulses  for  which  instantaneous  frequency  and  amplitude  become  ill-defined  quantities 
as  they  can  change  significantly  within  a  single  wave  cycle. 

Consider  the  time-domain  differential  form  of  the  Maxwell  equations 


d 

e(x,y,z)—E  +  a(x^y^z)E  =  vxh-j, 

d 

n(x,y,z)—n  + a*(x,y,z)n  =  -v  x  e, 


(28) 

(29) 


where  material  parameters  e,  /r,  a,  and  a*  are  the  electric  permittivity,  magnetic  permeability, 
electric  conductivity,  and  magnetic  conductivity,  respectively.  J  is  the  electric  current  density.  In 
inhomogeneous  media,  they  are  functions  of  spatial  variables  x ,  y,  and  z. 

There  are  two  ways  to  represent  the  Maxwell  equations  in  moving  frame:  the  Eulerian  and  the 
Lagrangian  approaches.  In  Eulerian  frame,  the  Maxwell  equations  stay  the  same  while  the  material 
parameters  are  functions  of  space  and  time: 


d 

e{x,y,z,t)—E  +  a(x,y,z,t)E  =  VxH-J, 

d 

V(x,  y,  z,  t)— H  +  a*  (x,  y,  z,  t) H  =  -V  x  E. 


(30) 

(31) 


In  Lagrangian  frame,  a  transformation  x'  =  x  —  ct  is  applied  where  x  is  assumed  to  be  the  direction 
of  moving  frame,  so  the  Maxwell  equations  have  additional  advection  terms: 


d  <9E 

e{x\y,z)-E-c—  +  cr(x',y,z)E  =  VxH-J, 
d  <9H 

^{x',y,z)—U-c—+a*(x',y,z)ll  =  -V  X  E. 


(32) 

(33) 
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Figure  10:  Simulation  of  optical  pulse  passing  through  a  dielectric  particle  using  stationary  and 
moving  FDTD  methods. 


We  have  studied  the  2D  moving  window  FDTD  in  Eulerian  frame  for  linear  media.  We  apply 
Berenger’s  split-field  approach  to  implement  the  PML  boundary  conditon.  Assume  the  moving 
window  travels  with  the  speed  of  light  in  x  direction  in  Eulerian  frame,  the  split-field  Maxwell 
equations  are  given  by 


d 

e(x,y,t)—Ex  +  ax(x,y,t)Ex 

dHz 

dy 

J Xl 

(34) 

d 

e(x>  ^oiEy  +  ay<yXl  y ’  ^Ey 

dHz 

dx 

—  Jyi 

(35) 

Q 

h( x,y,t)—Hzx  +  a*(x,  y,  t)Hzx 

Ctf  H 

1 

II 

(36) 

d 

^{x,y,t)—HZy  +  <Jy(x,  y,  t)Hzy 

II 

(37) 

HZx  Hzy 

=  Hz 

(38) 

For  Eulerian  approach,  the  Maxwell  equations  are  solved  by  the  standard  Yee  FDTD  method. 
The  difference  between  the  moving  and  the  stationary  FDTD  method  is  that  the  time  dependent 
material  parameters  need  to  be  evaluated  at  every  time  step.  Every  time  step,  the  moving  window 
is  moved  by  one  cell  toward  the  pulse  propagation  direction  if  the  center  of  the  main  pulse  has 
moved  out  of  the  center  cell,  otherwise  the  window  is  not  moved.  We  compare  our  moving  window 
FDTD  results  with  the  standard  FDTD  method  in  stationary  frame.  As  shown  in  Fig.  10,  good 
agreement  between  the  moving  and  stationary  FDTD  is  obtained  for  the  pulse  propagation  through 
a  dielectric  and  dispersive  cylinder.  The  error  is  less  than  1%  after  the  pulse  has  been  propagated 
for  a  relative  long  distance. 
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Figure  11:  Simulation  of  optical  pulse  passing  through  a  metal  (dispersive)  particle  using  stationary 
and  moving  FDTD  methods. 

6  Ground  Penetrating  Radar  simulations  by  Adaptive  Mesh  Re¬ 
finement  FDTD  method 

The  FDTD  method  has  been  applied  in  a  wide  range  of  applications  [3,  21,  22],  including  antenna, 
microwave  circuits,  geophysics,  optics,  etc.  The  Ground  Penetrating  Radar  (GPR)  is  a  popular  and 
efficient  nondestructive  electromagnetic  device  for  high-resolution  imaging  of  the  shallow  subsurface 
on  the  earth,  and  man-made  structures.  It  is  a  technique  that  has  been  employed  in  many  fields 
such  as  engineering,  geology,  and  environmental  studies.  For  example,  GPR  can  be  used  to  locate 
buried  utilities  and  to  measure  snow/ice  thickness.  It  can  also  be  used  as  wall  penetrating  radar 
to  ‘see’  through  the  wall.  GPR  systems  are  implemented  by  sending  a  pulse  into  a  material  via 
a  transmitter.  An  integrated  device  records  the  strength  and  time  required  for  the  return  of  any 
reflected  signals.  Subsurface  variations  will  create  reflections  that  are  picked  up  by  the  system 
and  recorded.  These  reflections  are  produced  by  a  variety  of  shapes  and  material,  which  in  turn 
are  used  to  form  images.  Under  favorable  conditions,  GPR  can  provide  very  accurate  information 
concerning  the  nature  of  various  subsurfaces  and  buried  objects. 

We  have  applied  the  Adaptive  Mesh  Refinement  (AMR)  FDTD  method  to  the  GPR  system. 
The  AMR  method  allows  the  simulation  of  very  small  object  with  fine  mesh  inside  a  large  region 
with  coarse  mesh.  Let  EV:„  and  E™  f  denote  the  field  values  defined  on  coarse  and  fine  meshes, 
respectively.  The  AMR  FDTD  algorithm  updating  the  Maxwell  equations  in  the  following  steps: 

1.  Update  E "+1  by  At  from  E™c. 

2.  Update  E'‘^  in  the  interior  of  the  fine  mesh  by  from  E™  j. 

n+i 

3.  Compute  E,  on  coarse/fine  boundary  through  interpolation  of  neighboring  coarse  mesh 
E™c  and  E™^1  electric  field  values. 

4.  Update  Hx  ^  and  H  ^  everywhere  on  the  fine  mesh  by  from  Hx  ^  and  Hy  j 7  respec- 
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Figure  12:  FDTD  simulation  of  GPR.  (a)  The  red  and  blue  circles  are  the  transmitter  and  receiver 
respectively,  (b)  A  pulse  is  emitted  at  multiple  locations  along  the  ground  and  the  data  from 
reflected  waves  from  underground  object  are  recorded,  (c):  FDTD  result  with  coarse  mesh  size  A 
everywhere,  (d):  FDTD  result  with  fine  mesh  size  A/2  everywhere,  (e):  AMR-FDTD  result  with 
fine  mesh  size  A/2  near  the  object  and  coarse  mesh  size  A  elsewhere. 


tively. 

n~\~  —  ti~\~  —  el~\~  —  ti~\~  — 

5.  Compute  HXjC2  and  HVtC2  through  interpolation  of  neighboring  fine  mesh  Htj4,  Hx  f 4 , 
n-|--  n+- 

H  f  4  and  H  ^  4  magnetic  field  values. 


6.  Update  F+ji1  in  the  interior  of  the  fine  mesh  by  +  from  u"^2. 

7.  Compute  77/ j1  on  coarse/fine  boundary  through  interpolation  of  neighboring  coarse  mesh 
17+ 1  electric  field  values. 

71+ -  71+ -  71+ -  71+ - 

8.  Update  Hx  +  and  Hy  +  everywhere  on  the  fine  mesh  from  Hx  +  and  Hy  ^ 4  respectively. 

9.  Compute  E through  interpolation  of  neighboring  fine  mesh  E™^1  field  values. 

10.  Update  Hx^2  and  Hy^2  by  At  from  Hxx2  and  Hy^2 . 


The  GPR  simulation  set  up  and  results  are  shown  in  Fig.  12.  We  compare  the  results  of  the 
standard  and  the  AMR  FDTD  methods.  For  AMR  simulation,  the  mesh  is  refined  in  a  small  region 
near  the  circular  object.  For  standard  FDTD  simulations,  the  fine  mesh  result  is  much  clearer  than 
that  of  the  coarse  mesh.  In  the  refined  region,  the  AMR  result  is  very  similar  to  that  of  the  fine 
mesh  FDTD  result,  but  the  simulation  takes  much  less  memory  and  CPU  time.  Therefore,  our 
AMR  FDTD  solver  provides  an  efficient  tool  in  GPR  simulation. 
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7  Conclusion  and  future  work 


The  main  objective  of  this  effort  is  the  development  of  stable,  accurate  and  efficient  Maxwell 
solvers.  The  proposed  activities  focus  on  mathematical  studies  of  the  key  unresolved  issues  in 
the  Finite-Difference  Time-Domain  (FDTD)  electromagnetic  simulations.  In  summary,  our  main 
achievements  and  discovers  are  listed  below. 

•  A  new  way  of  extending  the  subpixel  smoothing  FDTD  algorithm  to  dielectric/dispersive 
material  interface  was  invented  by  using  local  coordinate  rotation.  As  our  method  does 
require  split  fields,  better  efficiency  has  been  achieved  in  comparison  to  the  previous  work. 

•  A  novel  stable  anisotropic  overlapping  Yee  FDTD  Maxwell  solver  was  proposed  for  problem 
involving  complex  material  interface,  such  as  the  electromagnetic  invisibility  cloaking  devices. 
Previous  proposed  conventional  anisotropic  FDTD  method  suffers  from  the  late-time  insta¬ 
bility  problem  in  such  cases  due  to  the  extrapolation  near  the  material  interface.  In  contrast, 
our  method  relies  on  the  overlapping  cells  to  provide  the  collocated  field  values  without  any 
interpolation  or  extrapolation. 

•  A  locally  non-orthogonal  OY  FDTD  algorithm  was  developed  to  compute  optical  forces  on 
nanoparticles.  Together  with  diagonal  split  cell  mesh,  the  OY  approach  attains  second-order 
convergence  for  structures  with  smooth  curved  surface  and  with  sharp  corners. 

•  A  moving  window  FDTD  method  with  PML  boundaries  has  been  deveolped  to  model  high 
intensity  optical  pulse  propagation  over  long  distance  (millions  of  wavelength).  The  algorithm 
is  based  on  the  Eulerian  formulation.  Numerial  tests  show  that  the  moving  window  method 
is  accurate  in  comparison  to  the  stationary  FDTD  method. 

•  Adaptive  Mesh  Refinement  FDTD  method  has  been  applied  to  study  the  ground  penetrating 
radar  devices.  This  work  was  included  in  the  MS  thesis  of  a  graduate  student. 

In  the  future,  we  will  continue  our  investigation  on  new  algorithms  in  numerical  solutions  to 
Maxwell’s  equations  with  applications  to  linear  and  nonlinear  optics,  photonics,  EM  theory,  and 
metamaterials. 
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